rgrass7, raster, leaflet, leafem), fueron en el bloque anterior al ejecutarse el código contenido en el archivo extraer-red-de-drenaje-con-r-stream.Rmd. Igualmente, dicho bloque de código creó todos los objetos necesarios para realizar este tutorial.execGRASS(
'g.list',
flags = 't',
parameters = list(
type = c('raster', 'vector')
)
)
## raster/LfpNetwork-flds-Oz
## raster/LfpNetwork-outlet-Oz
## raster/MASK
## raster/accum-de-rwshed
## raster/aspect
## raster/basins
## raster/dem
## raster/drainage-dir-de-rstr
## raster/drainage-dir-de-rwshed
## raster/half-basins
## raster/order-hack-gravelius
## raster/order-horton
## raster/order-shreve
## raster/order-strahler
## raster/order-topology
## raster/ozama-basin
## raster/ozama-stream-de-rstr
## raster/pcurv
## raster/r-stream-basins-1
## raster/r-stream-basins-2
## raster/r-stream-basins-3
## raster/r-stream-basins-4
## raster/r-stream-basins-5
## raster/r-stream-basins-6
## raster/r-stream-basins-7
## raster/slope
## raster/stream-de-rwshed
## raster/tcurv
## vector/LfpNetwork_lfp_all_final_ozm
## vector/LfpNetwork_lfp_ozm
## vector/LfpNetwork_outlet_Oz
## vector/LfpNetwork_tributaries_ozm
## vector/LfpNetwork_tributaries_preconf_ozm
## vector/c_ozama
## vector/dem_extent
## vector/order_all
## vector/ozama_basin
## vector/ozama_stream_de_rstr
## vector/r_stream_basins_1
## vector/r_stream_basins_2
## vector/r_stream_basins_3
## vector/r_stream_basins_4
## vector/r_stream_basins_5
## vector/r_stream_basins_6
## vector/r_stream_basins_7
execGRASS(
"r.stream.extract",
flags = c('overwrite','quiet'),
parameters = list(
elevation = 'dem',
threshold = 80,
direction = 'drainage-dir-de-rstr'
)
)
execGRASS(
"r.stream.order",
flags = c('overwrite','quiet'),
parameters = list(
stream_rast = 'ozama-stream-de-rstr',
direction = 'drainage-dir-de-rstr',
elevation = 'dem',
accumulation = 'accum-de-rwshed',
stream_vect = 'order_all',
strahler = 'order-strahler',
horton = 'order-horton',
shreve = 'order-shreve',
hack = 'order-hack-gravelius',
topo = 'order-topology'
)
)
## Warning in execGRASS("r.stream.order", flags = c("overwrite", "quiet"), : The command:
## r.stream.order --overwrite --quiet stream_rast=ozama-stream-de-rstr direction=drainage-dir-de-rstr elevation=dem accumulation=accum-de-rwshed stream_vect=order_all strahler=order-strahler horton=order-horton shreve=order-shreve hack=order-hack-gravelius topo=order-topology
## produced at least one warning during execution:
## WARNING: Vector map <order_all> already exists and will be overwritten
## WARNING: Vector map <order_all> already exists and will be overwritten
execGRASS(
'g.list',
flags = 't',
parameters = list(
type = c('raster', 'vector')
)
)
## raster/LfpNetwork-flds-Oz
## raster/LfpNetwork-outlet-Oz
## raster/MASK
## raster/accum-de-rwshed
## raster/aspect
## raster/basins
## raster/dem
## raster/drainage-dir-de-rstr
## raster/drainage-dir-de-rwshed
## raster/half-basins
## raster/order-hack-gravelius
## raster/order-horton
## raster/order-shreve
## raster/order-strahler
## raster/order-topology
## raster/ozama-basin
## raster/ozama-stream-de-rstr
## raster/pcurv
## raster/r-stream-basins-1
## raster/r-stream-basins-2
## raster/r-stream-basins-3
## raster/r-stream-basins-4
## raster/r-stream-basins-5
## raster/r-stream-basins-6
## raster/r-stream-basins-7
## raster/slope
## raster/stream-de-rwshed
## raster/tcurv
## vector/LfpNetwork_lfp_all_final_ozm
## vector/LfpNetwork_lfp_ozm
## vector/LfpNetwork_outlet_Oz
## vector/LfpNetwork_tributaries_ozm
## vector/LfpNetwork_tributaries_preconf_ozm
## vector/c_ozama
## vector/dem_extent
## vector/order_all
## vector/ozama_basin
## vector/ozama_stream_de_rstr
## vector/r_stream_basins_1
## vector/r_stream_basins_2
## vector/r_stream_basins_3
## vector/r_stream_basins_4
## vector/r_stream_basins_5
## vector/r_stream_basins_6
## vector/r_stream_basins_7
order <- readVECT('order_all')
order4326 <- spTransform(order, CRSobj = CRS("+init=epsg:4326"))
leaflet() %>%
addProviderTiles(providers$Stamen.Terrain, group = 'terrain') %>%
addPolylines(
data = order4326, weight = 3, opacity = 0.7, group = 'order',
label = ~as.character(strahler),
highlightOptions = highlightOptions(color = "white",
weight = 5, bringToFront = F, opacity = 1),
labelOptions = labelOptions(noHide = F,
style = list(
"font-size" = "8px",
"background" = "rgba(255, 255, 255, 0.5)",
"background-clip" = "padding-box",
"padding" = "1px"))) %>%
leafem::addHomeButton(extent(order4326), 'Ver todo') %>%
addLayersControl(
overlayGroups = c('terrain','order'),
options = layersControlOptions(collapsed=FALSE))
orden_de_red <- leaflet() %>%
addProviderTiles(providers$Stamen.Terrain, group = 'terrain') %>%
addPolylines(
data = order4326, weight = order4326$strahler*1.5, opacity = 0.7, group = 'order',
label = ~as.character(strahler),
highlightOptions = highlightOptions(color = "white",
weight = 5, bringToFront = F, opacity = 1),
labelOptions = labelOptions(noHide = F)) %>%
leafem::addHomeButton(extent(order4326), 'Ver todo') %>%
addLayersControl(
overlayGroups = c('terrain','order'),
options = layersControlOptions(collapsed=FALSE))
orden_de_red
orden_de_red %>% mapview::mapshot(file = 'orden_de_red_salida.png')
#Estadísticas para obtener los valores mínimo y máximo del orden de red de Strahler
rinfo.ordstra <- execGRASS(
'r.info',
flags = 'r',
parameters = list(
map = 'order-strahler'
)
)
## min=1
## max=7
#Órdenes de red mínimo y máximo
minmaxord <- as.numeric(
stringr::str_extract_all(
attributes(rinfo.ordstra)$resOut,
"[0-9]+"
)
)
minmaxord
## [1] 1 7
sapply(
min(minmaxord):max(minmaxord),
function(x){
execGRASS(
"r.stream.basins",
flags = c('overwrite','c','quiet'),
parameters = list(
direction = 'drainage-dir-de-rstr',
stream_rast = 'order-strahler',
cats = as.character(x),
basins = paste0('r-stream-basins-',x)
)
)
execGRASS(
"r.to.vect",
flags=c('overwrite','quiet'),
parameters = list(
input = paste0('r-stream-basins-',x),
output = paste0('r_stream_basins_',x),
type = 'area'
)
)
}
)
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-1 output=r_stream_basins_1 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_1> already exists and will be
## overwritten
## WARNING: Vector map <r_stream_basins_1> already exists and will be
## overwritten
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-2 output=r_stream_basins_2 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_2> already exists and will be
## overwritten
## WARNING: Vector map <r_stream_basins_2> already exists and will be
## overwritten
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-3 output=r_stream_basins_3 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_3> already exists and will be
## overwritten
## WARNING: Vector map <r_stream_basins_3> already exists and will be
## overwritten
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-4 output=r_stream_basins_4 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_4> already exists and will be
## overwritten
## WARNING: Vector map <r_stream_basins_4> already exists and will be
## overwritten
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-5 output=r_stream_basins_5 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_5> already exists and will be
## overwritten
## WARNING: Vector map <r_stream_basins_5> already exists and will be
## overwritten
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-6 output=r_stream_basins_6 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_6> already exists and will be
## overwritten
## WARNING: Vector map <r_stream_basins_6> already exists and will be
## overwritten
## Warning in execGRASS("r.to.vect", flags = c("overwrite", "quiet"), parameters = list(input = paste0("r-stream-basins-", : The command:
## r.to.vect --overwrite --quiet input=r-stream-basins-7 output=r_stream_basins_7 type=area
## produced at least one warning during execution:
## WARNING: Vector map <r_stream_basins_7> already exists and will be
## overwritten
## WARNING: Vector map <r_stream_basins_7> already exists and will be
## overwritten
## [1] 0 0 0 0 0 0 0
sapply(
min(minmaxord):max(minmaxord),
function(x){
assign(
paste0('orden', x),
spTransform(readVECT(paste0('r_stream_basins_',x),driver = 'SQLite'), CRSobj = CRS("+init=epsg:4326")),
envir = .GlobalEnv)
}
)
paleta <- RColorBrewer::brewer.pal(12, 'Set3')
cuencas_y_orden_de_red <- leaflet() %>%
addProviderTiles(providers$Stamen.Terrain, group = 'terrain') %>%
addPolygons(data = orden7, stroke = T, weight = 2,
color = ~paleta, fillOpacity = 0.4, group = 'O7') %>%
addPolygons(data = orden6, stroke = T, weight = 2,
color = ~paleta, fillOpacity = 0.4, group = 'O6') %>%
addPolygons(data = orden5, stroke = T, weight = 2,
color = ~paleta, fillOpacity = 0.4, group = 'O5') %>%
addPolygons(data = orden4, stroke = T, weight = 2,
color = ~paleta, fillOpacity = 0.4, group = 'O4') %>%
addPolygons(data = orden3, stroke = T, weight = 2,
color = ~paleta, fillOpacity = 0.4, group = 'O3') %>%
addPolygons(data = orden2, stroke = T, weight = 2,
color = ~paleta, fillOpacity = 0.4, group = 'O2') %>%
addPolygons(data = orden1, stroke = T, weight = 2,
color = ~paleta, fillOpacity = 0.4, group = 'O1') %>%
addPolylines(
data = order4326, weight = order4326$strahler*1.5,
opacity = 0.7, group = 'str_order') %>%
leafem::addHomeButton(extent(order4326), 'Ver todo') %>%
addLayersControl(
overlayGroups = c('terrain','O1','O2','O3','O4','05','06','07','str_order'),
options = layersControlOptions(collapsed=FALSE))
cuencas_y_orden_de_red
cuencas_y_orden_de_red %>% mapview::mapshot(file = 'cuencas_y_orden_de_red_salida.png')
execGRASS(
"r.stream.stats",
flags = c('overwrite','quiet','o'),
parameters = list(
stream_rast = 'order-strahler',
direction = 'drainage-dir-de-rstr',
elevation = 'dem',
output = 'ozama_stats.txt'
)
)
file.show('ozama_stats.txt')
d <- read.csv("ozama_stats.txt", skip=1, header=TRUE)
d
## order num_of_streams avg_length avg_area avg_slope avg_grad
## 1 1 1328 1.157222 1.531833 0.019501 0.013230
## 2 2 274 2.787560 6.990865 0.013704 0.007789
## 3 3 63 8.174857 36.891303 0.009048 0.004075
## 4 4 13 18.228154 180.416363 0.007264 0.001903
## 5 5 5 21.793405 494.576211 0.004899 0.000828
## 6 6 2 17.109400 1364.385565 0.002296 0.000297
## 7 7 1 23.929262 3376.766977 0.005338 0.000190
## avg_elev.diff sum_length sum_area
## 1 15.912760 1536.79041 2034.274
## 2 21.873087 763.79156 1915.497
## 3 39.441891 515.01597 2324.152
## 4 45.280315 236.96600 2345.413
## 5 12.680176 108.96702 2472.881
## 6 5.009900 34.21880 2728.771
## 7 4.555471 23.92926 3376.767
plot(num_of_streams~order, data=d, log="y")
mod <- lm(log10(num_of_streams)~order, data=d)
abline(mod)
text(2, 20, 'logN=2.064-0.544u')
rb <- 1/10^mod$coefficients[[2]]
rb
## [1] 3.361632
execGRASS(
"r.stream.stats",
flags = c('overwrite','quiet'),
parameters = list(
stream_rast = 'order-strahler',
direction = 'drainage-dir-de-rstr',
elevation = 'dem',
output = 'ozama_stats_expanded.txt'
)
)
file.show('ozama_stats_expanded.txt')
rbp <- mean(d$num_of_streams[-length(d$num_of_streams)]/d$num_of_streams[-1])
rbp
## [1] 3.523679
unlink_.gislock()